library(foreign)
library(coda)

setwd("E:/finalidealreplication")

## Load MCMC chains from one-cutpoint model
load("onecutpointdraws.RData")

## Load signature matrix
mpinc98.v4 <- read.csv("mpinc98.v4.csv", row.names=1)
sigs.98 <- apply(mpinc98.v4, 1, sum) 

## Load MP characteristics
mpchar <- read.csv("mpcharacteristics.csv", row.names=1)
mp98party <- mpchar$mp98party
mp98party.color <- character(length(mp98party))
mp98party.color <- "black"
mp98party.color[mp98party=="Lab"] <- "red"
mp98party.color[mp98party=="Con"] <- "blue"
mp98party.color[mp98party=="LDem"] <- "yellow"

## Calculate posterior means
mean.idealpoints <- apply(idealpoints,2,mean)
mean.statusquo <- apply(statusquo,2,mean)
mean.motion <- apply(motion,2,mean)
mean.propensity <- apply(propensity,2,mean)

## Reorder MPs by posterior mean ideal point estimates
order.idealpoints <- order(mean.idealpoints, decreasing = TRUE)


## Load vector of sponsor row IDs
sponsorchar <- read.csv("sponsors.csv", row.names=1)
sponsor.vector98 <- sponsorchar[,1]
sponsor.party <- numeric()
sponsor.ideal <- numeric()
for (i in 1:length(sponsor.vector98)){
sponsor.party[i] <- mp98party[sponsor.vector98[i]]
sponsor.ideal[i] <- mean.idealpoints[sponsor.vector98[i]]
}

sponsor.party <- as.factor(sponsor.party)
levels(sponsor.party) <- c("Con", "DU", "LDem", "Lab", "PC", "SDLP", "SNP", "UU", "WW")





## Construct Table 1 of Kellermann (2012)
library(MASS)
summary(rlm((mpchar[mp98party=="Lab",5]) ~ mean.idealpoints[mp98party=="Lab"], method="MM"))
summary(rlm((mpchar[which(mp98party=="LDem"),5]) ~ mean.idealpoints[which(mp98party=="LDem")], method="MM"))
summary(rlm((mpchar[mp98party=="Con",5]) ~ mean.idealpoints[mp98party=="Con"], method="MM"))



## Construct Figure 1 of Kellermann (2012)
pdf("idhists.pdf", width=9, height=3)
plot(density(idealpoints[,which(rownames(mpinc98.v4)=="Radice/Giles")]), xlim=c(-1.5, 1), ylim=c(0,10), main="", xlab="Estimated ideal point")
lines(density(idealpoints[,which(rownames(mpinc98.v4)=="Abbott/Diane")]))
lines(density(idealpoints[,which(rownames(mpinc98.v4)=="Livingstone/Ken")]))
lines(density(idealpoints[,which(rownames(mpinc98.v4)=="Bell/Martin")]))
lines(density(idealpoints[,which(rownames(mpinc98.v4)=="Bottomley/Virginia")]))
lines(density(idealpoints[,which(rownames(mpinc98.v4)=="Chope/Christopher")]))

text(.45, 5, "Christopher Chope")
text(0, 6.5, "Virginia Bottomley")

text(-1.1, 9, "Ken Livingstone")
text(-.85, 7.5, "Diane Abbott")

text(-.6, 5, "Giles Radice")
text(-.35, 8, "Martin Bell")
dev.off()


## Construct Figure 2 of Kellermann (2012)
pdf("idealcost2.pdf")
plot(mean.idealpoints[which(mp98party=="Lab")],-1*mean.propensity[which(mp98party=="Lab")], pch=1, col="black", xlab = "Estimated ideal points", ylab="Estmated cost of signing", xlim=range(mean.idealpoints), ylim=c(0,max(-1*mean.propensity)) )
points(mean.idealpoints[which(mp98party=="Con")],-1*mean.propensity[which(mp98party=="Con")], pch=2, col="black")
points(mean.idealpoints[which(mp98party=="LDem")],-1*mean.propensity[which(mp98party=="LDem")], pch=7, col="black")
##identify(mean.idealpoints, mean.propensity, labels=rownames(mpinc98.v4))
legend(-1.5, 0.75, c("Labour", "Conservative", "LibDem"), pch=c(1,2,7), bty="n")
dev.off()




## Create smaller sample of posterior draws for computational purposes
short.idealpoints <- idealpoints[4*1:4000, ]

## Construct posterior ranks
idealpoints.rank <- matrix(NA, nrow(short.idealpoints), ncol(short.idealpoints))
for(i in 1:4000){
idealpoints.rank[i,] <- rank(short.idealpoints[i,])
}


## Construct Figure 3 of Kellermann (2012)
pdf("idealinterval2.pdf", height=9)
par(mfrow=c(3,2), oma=c(0,2,2,0), mar=c(4,4,2,1)+.1)
plot(1, 1, xlim=c(-1.5,1.5), ylim=c(0,550), type="n", xlab="95% credible intervals", ylab="Rank")
for(i in which(mp98party.color[order.idealpoints]=="red")){
segments( quantile(short.idealpoints[,order.idealpoints[i]], .025), i, quantile(short.idealpoints[,order.idealpoints[i]], .975),i, col="darkgrey")
}

plot(1, 1, xlim=c(0,550), ylim=c(0,550), type="n", xlab="95% credible intervals", ylab="Rank")
for(i in which(mp98party.color[order.idealpoints]=="red")){
segments( quantile(idealpoints.rank[,order.idealpoints[i]], .025), i, quantile(idealpoints.rank[,order.idealpoints[i]], .975),i, col="darkgrey")
}

plot(1, 1, xlim=c(-1.5,1.5), ylim=c(0,550), type="n", xlab="95% credible intervals", ylab="Rank")
for(i in which(mp98party.color[order.idealpoints]=="yellow")){
segments( quantile(short.idealpoints[,order.idealpoints[i]], .025), i, quantile(short.idealpoints[,order.idealpoints[i]], .975),i, col="darkgrey")
}

plot(1, 1, xlim=c(0,550), ylim=c(0,550), type="n", xlab="95% credible intervals", ylab="Rank")
for(i in which(mp98party.color[order.idealpoints]=="yellow")){
segments( quantile(idealpoints.rank[,order.idealpoints[i]], .025), i, quantile(idealpoints.rank[,order.idealpoints[i]], .975),i, col="darkgrey")
}

plot(1, 1, xlim=c(-1.5,1.5), ylim=c(0,550), type="n", xlab="95% credible intervals", ylab="Rank")
for(i in which(mp98party.color[order.idealpoints]=="blue")){
segments( quantile(short.idealpoints[,order.idealpoints[i]], .025), i, quantile(short.idealpoints[,order.idealpoints[i]], .975),i, col="darkgrey")
}

plot(1, 1, xlim=c(0,550), ylim=c(0,550), type="n", xlab="95% credible intervals", ylab="Rank")
for(i in which(mp98party.color[order.idealpoints]=="blue")){
segments( quantile(idealpoints.rank[,order.idealpoints[i]], .025), i, quantile(idealpoints.rank[,order.idealpoints[i]], .975),i, col="darkgrey")
}

mtext("Estimated ideal points", side=3, at=.28,outer=T, font=2)
mtext("Estimated rank", side=3, at=.77,outer=T, font=2)
mtext("Conservative", side=2, at=.18,outer=T, cex=.85)
mtext("Liberal Democrat", side=2, at=.51,outer=T, cex=.85)
mtext("Labour", side=2, at=.85,outer=T, cex=.85)

dev.off()


################

load("twocutpointdraws.RData")

## Calculate posterior means from two-cutpoint model
mean.ideal.sal <- apply(ideal.sal,2,mean)
mean.salience.sal <- apply(salience.sal,2,mean)
mean.motion.sal <- apply(motion.sal,2,mean)
mean.propensity.sal <- apply(propensity.sal,2,mean)


## Repeats robust linear regression using estimates from two-cutpoint model
summary(rlm((mpchar[mp98party=="Lab",5]) ~ mean.ideal.sal[mp98party=="Lab"], method="MM"))
summary(rlm((mpchar[which(mp98party=="LDem"),5]) ~ mean.ideal.sal[which(mp98party=="LDem")], method="MM"))
summary(rlm((mpchar[mp98party=="Con",5]) ~ mean.ideal.sal[mp98party=="Con"], method="MM"))

## Constructs Figure 1 of supplemental information for Kellermann (2012)
pdf("E:/unfold/suppgraph3.pdf", width=8, height=4)
par(mfrow=c(1,2))
plot(mean.idealpoints, sigs.98, xlab="Estimated ideal point", ylab="Number of signatures" )
title("(a) One-cutpoint model")

plot(mean.ideal.sal, sigs.98, xlab="Estimated ideal point", ylab="Number of signatures")
title("(b) Two-cutpoint model")
dev.off()


## Construct Figure 2 of supplemental information for Kellermann (2012)
pdf("E:/unfold/suppgraph2.pdf", width=8, height=4)
par(mfrow=c(1,2))

plot(mean.idealpoints[mp98party=="Lab"], mpchar[mp98party=="Lab",5], col="darkgrey", xlab="Estimated ideal points", ylab="Labor voting support", xlim=c(-2, 1.0), ylim=c(0.8,1))
lines(loess.smooth( mean.idealpoints[mp98party=="Lab"], mpchar[mp98party=="Lab",5]), lwd=2)

#points(mean.idealpoints[which(mp98party=="LDem")], mpchar[which(mp98party=="LDem"),5], col="darkgrey", pch=7)
#lines(loess.smooth( mean.idealpoints[which(mp98party=="LDem")], mpchar[which(mp98party=="LDem"),5]), lwd=2)

#points(mean.idealpoints[mp98party=="Con"], mpchar[mp98party=="Con",5], col="darkgrey", pch=2)
#lines(loess.smooth( mean.idealpoints[mp98party=="Con"], mpchar[mp98party=="Con",5]), lwd=2)

title("(a) One-cutpoint model")

plot(mean.ideal.sal[mp98party=="Lab"], mpchar[mp98party=="Lab",5], col="darkgrey", xlab="Estimated ideal points", ylab="Labor voting support", xlim=c(-3, 2.0), ylim=c(0.8,1))
lines(loess.smooth( mean.ideal.sal[mp98party=="Lab"], mpchar[mp98party=="Lab",5]), lwd=2)

#points(mean.ideal.sal[which(mp98party=="LDem")], mpchar[which(mp98party=="LDem"),5], col="darkgrey", pch=7)
#lines(loess.smooth( mean.ideal.sal[which(mp98party=="LDem")], mpchar[which(mp98party=="LDem"),5]), lwd=2)

#points(mean.ideal.sal[mp98party=="Con"], mpchar[mp98party=="Con",5], col="darkgrey", pch=2)
#lines(loess.smooth( mean.ideal.sal[mp98party=="Con"], mpchar[mp98party=="Con",5]), lwd=2)

title("(b) Two-cutpoint model")
dev.off()


## Constructs Figure 3 of supplemental information for Kellermann (2012)
pdf("motionshift.pdf", width=9, height=3)

par(mfrow=c(1, 3))
plot(mean.motion[which(sponsor.party=="Lab")],mean.statusquo[which(sponsor.party=="Lab")], pch=1, col="red", xlab="Estimated motion positions", ylab = "Estimated status quo positions", ylim=c(-2.5, 2.5), xlim=c(-2.5, 2.5), ##main = "Figure 3: Posterior means of motion parameters"
, type="n")
arrows(mean.statusquo[which(sponsor.party=="Lab")], mean.statusquo[which(sponsor.party=="Lab")], mean.motion[which(sponsor.party=="Lab")],mean.statusquo[which(sponsor.party=="Lab")], length=.1, col="darkgrey")
title("Labour-sponsored EDMs")


plot(mean.motion[which(sponsor.party=="Lab")],mean.statusquo[which(sponsor.party=="Lab")], pch=1, col="red", xlab="Estimated motion positions", ylab = "Estimated status quo positions", ylim=c(-2.5, 2.5), xlim=c(-2.5, 2.5), ##main = "Figure 3: Posterior means of motion parameters"
, type="n")
arrows(mean.statusquo[which(sponsor.party=="LDem")], mean.statusquo[which(sponsor.party=="LDem")], mean.motion[which(sponsor.party=="LDem")],mean.statusquo[which(sponsor.party=="LDem")], length=.1, col="darkgrey")
title("Liberal Democrat-sponsored EDMs")

plot(mean.motion[which(sponsor.party=="Lab")],mean.statusquo[which(sponsor.party=="Lab")], pch=1, col="red", xlab="Estimated motion positions", ylab = "Estimated status quo positions", ylim=c(-2.5, 2.5), xlim=c(-2.5, 2.5), ##main = "Figure 3: Posterior means of motion parameters"
, type="n")
arrows(mean.statusquo[which(sponsor.party=="Con")], mean.statusquo[which(sponsor.party=="Con")], mean.motion[which(sponsor.party=="Con")],mean.statusquo[which(sponsor.party=="Con")], length=.1, col="darkgrey")
title("Conservative-sponsored EDMs")
dev.off()

## Constructs Figure 4 of supplemental information for Kellermann (2012)
pdf("costprecision.pdf")
plot(-mean.propensity, apply(idealpoints, 2, sd), xlab="Estimated cost of signing", ylab="Standard devation of estimated ideal point")
dev.off()


## Construct table of idealpoint estimates (Table 1 of supplemental information for Kellermann (2012)

neworder <- order.idealpoints

quantile.idealpoints <- apply(idealpoints, 2, quantile, c(.025, .975))

table.idealpoints <- cbind(rownames(mpinc98.v4)[neworder], as.character(mp98party[neworder]),round(mean.idealpoints[neworder], 2), round(t(quantile.idealpoints)[neworder,], 2))
rownames(table.idealpoints) <- NULL

write.table(table.idealpoints, sep="&", eol="\\\\\n", quote=F)

